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An optimized phonon approach for the numerical diagonalization of interacting electron-phonon 
systems is proposed. The variational method is based on an expansion in coherent states that leads 
to a dramatic truncation in the phonon space. The reliability of the approach is demonstrated 
for the extended Holstein model showing that different types of lattice distortions are present at 
intermediate electron-phonon couplings as observed in strongly correlated systems. The connection 
with the density matrix renormalization group is discussed. 
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The interaction of electrons with local lattice deforma- 
tions plays an important role in many different materi- 
als leading to a number of effects that range from colos- 
sal magnetoresistance in manganites0 to superconduc- 
tivity in fullerenes,0 from pseudo-gap in cupratesQ to 
Peierls instability in quasi-one-dimensional materials, 
only to mention a few. Although the theoretical study 
of the electron-phonon (e-ph) interaction has attracted 
constantly the interest of the scientific community, sev- 
eral aspects are still not fully understood and challeng- 
ing. Actually, analytical solutions of the e-ph problem are 
available only in weak and strong coupling asymptotic 
regimes, even for the simplest molecular crystal model 
introduced by Holstein. [f| Some insight in the problem 
for intermediate coupling comes from the dynamical field 
approach^ that provides exact results in the infinite di- 
mension limit and variational methods. In order to 
achieve a complete understanding of the e-ph interaction, 
numerical techniques as exact diagonalization and Quan- 
tum Monte Carlo have been exploited. In both cases the 
specific problem posed by the e-ph interaction is related 
to the presence of phonons that require an infinite dimen- 
sional Hilbert space. As a consequence, for instance, the 
direct application of Lanczos exact diagonalization has 
been limited to very small lattices and needs a somehow 
arbitrary truncation in the phonon number. In this con- 
text the introduction of an optimized phonon approach 
based on the analysis of the density matrixj^J has pro- 
duced an important improvement. The idea beyond this 
approach is to take advantage of the knowledge of the 
largest eigenvalues and eigenvectors of the site density 
matrix to select the phonon linear combinations that, at 
the best, can describe the system: the so-called opti- 
mized phonon basis (OPB). Unfortunately, the density 
matrix of the target states is not available a priori: it 
must be calculated in a self-consistent way together with 
the OPB. To this aim different strategies have been dis- 
cussed. In this report we wish to introduce a variational 
technique, based on an expansion in coherent states, that 
very much simplifies the selection of an optimized phonon 
basis and that does not require any truncation in the 



number of phonons. The method is able to provide ac- 
curate results for any e-ph coupling regime and for any 
value of the adiabatic ratio. As we will show in the case 
of the Holstein model, the proposed expansion provides 
states surprisingly close to the eigenvectors of the site 
density matrix corresponding to the highest probabilities 
and allows us to clarify some questions of the e-ph in- 
teraction that have been debated recently. In particular, 
the role of quantum lattice fluctuations in the Holstein 
model at half-filling and in the clustered phase is dis- 
cussed emphasizing the coexistence of very different lat- 
tice deformations that all contribute to the ground state. 

The model and the proposed optimized phonon basis. 
In the framework of e-ph systems the Holstein model has 
been very much studied. Here we discuss a simple gener- 
alization that includes nearest neighbor e-ph interaction 
terms. Furthermore, we consider spinless fermions in or- 
der to simulate at a very rough level0 a strong on-site 
electron repulsion. In the ID case the Hamiltonian is 
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where cj (cj) is the site electron creation (annihilation) 
operator and b\ creates (annihilates) a phonon in the 
site i. The phonon frequency is u>, the electron hopping 
between nearest neighbors is controlled by t while g and 
ge describe the strength of the on-site and nearest neigh- 
bor e-ph interactions, respectively. We assume K = 1. 

The "natural" basis in which we can describe Hamil- 
tonian Q is given by 
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where i run over the lattice sites, Vi is the electron state 
label that, for spinless electrons, can assume only two 
values and [ii labels the infinite phonon states on-site 
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i. Our aim is to achieve a very satisfying estimation of 
the ground state of fl} by using only a finite number of 
states. The first step is to replace the phonon states |/^) 
at site i with coherent states (CS) 
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where |0)^ p '* > is the phonon vacuum at the site i. This 
substitution does not introduce any truncation in the 
Hilbert space since, varying h in the complex plane, the 
local basis (0 is over-complete. The second step is to 
choose a finite number M of CS and the corresponding 
values h a (a = 1, M) that we have supposed not to 
depend on the site. The most efficient way to proceed 
is to fix the number of CS arbitrarily and choose the 
h a variationally calculating the ground state energy and 
minimizing it. In this way we obtain the "best" sub- 
space where to find the ground state. The method has 
the advantage to allow a systematic improvement of the 
approximation adding more and more CS and provides a 
quantitative estimation of the error introduced by a spe- 
cific truncation. Of course, the number of CS per site 
cannot grow as we wish since the size of the matrix to be 
diagonalized can become very large. For instance, if we 
choose two CS per site, the phonon matrix at fixed elec- 
tron configuration will be of size L — 2 N , where N is the 
number of the lattice sites. As we will show, two states 
per site are a good choice for a very accurate estimation 
of the ground state. 

The proposed optimized phonon basis can further be 
improved introducing a dependence on the electron state 
at site i since, on physical ground, the lattice deformation 
of site i is expected to depend very much on whether 
or not the site is occupied by an electron. This can be 
done acting with the operator S of the usual Lang-Firsov 
transformation^^ on each CS belonging to the chosen 
basis : 



S = e 9l(ni)+f(ni-(ni))](.bi-b\) 
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where / depends neither on site i nor on the specific CS 
to which S is applied. The value of / can be estimated 
variationally together with the set {h a }. Technically, it 
is easier first to work out the action of S transforming 
H to H' = S^HS and, then, to study the properties of 
H' in the subspace spanned by both the electron and the 
optimized phonon basis: 
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where h ai € {h a } and a, = f, ...,M. 

Single polaron features and site density matrix. In or- 
der to clarify the relation between the proposed approach 
and the optimized phonon states obtained within the 
density matrix method, we study a two-site one-electron 
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FIG. 1: The probability P of the N phonon state for the 
eigenvector corresponding to the largest eigenvalue of the site 
density matrix (/3 = f, occupied site). The system is made of 
two lattice sites and one electron. The numerically exact re- 
sults (stars) are compared to the results obtained within our 
approach for two (squares) and three (circles) CS. The values 
of the parameters are chosen in the intermediate coupling re- 
gion ("worst" case). In the inset: the second most important 
eigenvalue. 



system. In particular, we can construct the site den- 
sity matrixjg calculating its eigenvalues and eigenvec- 
tors. Both these quantities can be evaluated both within 
our approach and by an exact diagonalization method in 
the "natural" phonon basis. The comparison between the 
two results is contained in Fig.l where the two eigenvec- 
tors corresponding to the largest eigenvalues are plotted 
as a function of the phonon number for (3 — 1 (occupied 
site). It is clear that already with two CS the agreement 
is surprisingly good while with three CS the differences 
with the numerically exact results are completely negli- 
gible. The agreement is due to the fact that the exact 
site density matrix contains only a few relevant eigenval- 
ues and that the corresponding eigenvectors can be ap- 
proximated as linear combinations of CS. It is worthwhile 
noting that, roughly speaking, in the strong and interme- 
diate coupling regime, the two most important eigenvec- 
tors correspond to strong (weak) lattice deformations and 
weak (strong) lattice deformations, respectively. The two 
contributions are both important in the so-called inter- 
mediate regime while in the weak and strong regimes only 
one type of lattice deformation is relevant. Also for (3 = 
(empty site) the site deformation is built up through both 
weak and/or strong contributions; however, the method 
proposed is so flexible to allow different distortions for 
empty and occupied sites. These observations make clear 
why the generalized Lang Firsov method[ll|. based on a 
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TABLE I: The ground state energy E/t of a polaron for dif- 
ferent regimes by using 2, 3, 4 CS p er site is compared with re- 
sults provided by Wellein et al.jlSf. The system size is L — 10. 



g 2 ui/t E/t[15] 



2CS 



3CS 



4CS 



0.5 0.4 -2.06130 -2.06128 -2.06130 -2.06130 
5.0 0.4 -2.72827 -2.71365 -2.72777 -2.72828 
1.0 4.0 -5.02985 -5.02883 -5.02984 -5.02985 
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L=16 

i^p present work, tj/t=1.0 

A Ref. 112], w/t=1.0 

present work, u/t=0.5 

if present work, a/t=2.0 
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FIG. 2: The ground-state energy of a polaron as a function 
of the coupling constant A = g 2 u>/2t for different adiabatic 
ratios ui/t. The case co/t = 1 is compared with the data of 
Ref.0 



single CS, is not able to give satisfactory results in the in- 
termediate regime (Fig.l). One of the advantages of the 
proposed method is the possibility to drastically reduce 
the size of the matrix that represents the system Hamil- 
tonian (JTJ. We were able with little computational effort 
to calculate the ground state energy of a single polaron 
in a 16 site lattice by using two CS. The comparison with 
current estimations in the literature 0, 0] is excellent 
(Fig. 2 and Table 1) showing that two CS are enough to 
give a correct description in any regime. 

The half-filling case: anti- adiabatic regime and quan- 
tum fluctuations. It is well known that the ground state 
of the ID Holstein model exhibits a phase transition from 
a Luttinger liquid to a charge ordering state at half- 
filling, |13| but the role of the quantum fluctuations has 
not been fully explored. First of all we have checked that 
the total energy and the static structure factor calculated 
with our method are in excellent agreement with previ- 
ous estimations based on exact diagonalization obtained 
by Lanczos method. |l4| In particular, in the intermediate 
regime (the worst case) , our estimations based on two CS 
is 0.1% higher than the best estimation available. 15] In 



order to understand the role of the quantum fluctuations 
in the anti-adiabatic regime (lo 3> t) where quantum lat- 
tice effects are expected to play an important role, we 
have calculated P(X;x), the lattice displacement proba- 
bility distribution function (LDPDF), at a fixed electron 
configuration \. In Fig. 3 we show the results of the 
calculated LDPDF, relative to a specific site, associated 
to the most important electron configurations close to 
the boundary of the CO phase for lo = 5t and g 2 = 3.6. 
The full line gives the LDPDF of an occupied site in 
the CO electron configuration where alternating sites are 
occupied. As expected it is peaked at a large value of 
the lattice displacement (roughly X y/2Muj ~ 2g ~ 3.8, 
M being the ionic mass) since we are in the strong cou- 
pling regime. However, in Fig. 3 it is shown that there 
are less important electron configurations among which 
some contribute to the LDPDF with a peak at around 
Xq\/2M(jj ~ typical of empty sites. These less impor- 
tant electron configurations can be viewed as "defects" 
with respect to the "ideal" CO configuration. These 
defects are associated to very different lattice displace- 
ments compared to the ones characteristic of the dom- 
inant CO configuration, but are still important in the 
anti-adiabatic regime. We expect that this behavior is 
still present in the thermodynamic limit since the elec- 
tron configurations shown in Fig. 3 are relevant also in 
the broken symmetry phase of the thermodynamic limit. 
We note that the strong fluctuations represent the quan- 
tum counterpart of the thermal fluctuations found in the 
opposite adiabatic limit within the dynamical mean field 
approach (DMF).^|| Furthermore, we stress that, in the 
adiabatic limit, the quantum fluctuations behave in a 
very different way giving rise to "defects" whose lattice 
displacement is of the same type of the dominant CO 
configuration. For this reason, as shown in DMF ap- 
proach, the ground state LDPDF is characterized by a 
single peak. 

Evidence of quantum fluctuations in clustered phases 
at intermediate coupling. Recently the formation of clus- 
ters due to e-ph interaction has been discussed in con- 
nection to experimental evidences of charge ordering in 
cuprates and maneranites. [TR lisj One of the simplest 
models that supports cluster formation is the modified 
Holstein model of cq.Q. In the weak coupling regime 
the system is expected to form a Luttinger liquid, while, 
for e / the effective electron interaction, that is at- 
tractive for nearest neighbors, is responsible for the for- 
mation of electron clusters at strong coupling. As we 
go from strong to weak coupling, the cluster tends to 
break up and, in the intermediate coupling regime and 
for phonon and electron energy scales not well separated 
(uj ~ t), all the electron configurations contribute to the 
ground state without a dominant configuration (Fig. 4). 
In this regime, that is the most interesting both from 
the experimental and the theoretical point of view, our 
calculation shows that different electron configurations 
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FIG. 3: The LDPDF of the site indicated by the arrow for 
the most important electron configurations are plotted as a 
function of the lattice distortion in the anti-adiabatic region 
and in the charge-ordering regime. There are important elec- 
tron configurations characterized by very different distortions. 
The lattice distortions are measured in units of (2Mu)~ 1 ^ 2 . 



FIG. 4: The LDPDF of the site indicated by the arrow for 
two typical electron configurations are plotted as a function of 
the lattice distortion in the intermediate coupling regime for 
a system with nearest-neighbor e-ph interaction. The lattice 
distortions are measured in units of (2Mw) -1 ' 2 . The system 
is made of three electrons on L = f f lattice sites. 



can be associated to very different lattice deformations. 
This result stems out clearly from Fig. 4 where we have 
plotted the LDPDF of two electron configurations for 
A = 1.2 and e = 0.1. The site occupied by the cen- 
tral electron in the cluster configuration (a) is character- 
ized by a LDPDF peaked at high values of deformation 
while the site occupied by the central electron in config- 
uration (b) presents much weaker deformations typical 
of the weak coupling regime with larger spreading. The 
balance between effective hopping, reduced by the cluster 
formation, and e-ph energy gain can allow a larger lattice 
deformation energy. The presence of electron configura- 
tions with very different lattice deformations that con- 
tribute to the ground state almost at same level suggests 
the existence of significant quantum dynamical fluctua- 
tions triggered by e-ph interactions. We think that our 
result, although limited to a simple ID e-ph model, pro- 
vides support to a number of experimental evidences in 
cuprates [lJll2J| and manganites 1] where lattice distor- 
tions have been observed by infrared, neutron scattering 
and X-ray spectroscopy. 

Conclusions. In this report we have proposed a varia- 
tional approach, based on CS expansion, that is able to 
identify optimized phonon basis useful for studying e-ph 
models. The reliability of the model has been demon- 
strated for the extended Holstein model emphasizing the 
role of quantum fluctuations. In particular, it has been 
pointed out that very different lattice distortions coexist 
at intermediate e-ph coupling. 
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